%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fourier grid for solving a LiRb potential and mean distance.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function EstadoLivre



Constant;

	figure(1);						clf;

	cnst=e0^2/4/pi/E0/aB/2/pi/hP/100/c0;											%conversion a.u. to cm^-1





lvl='(1)1Pi';		eval(['[Rp,Vp]=LiRbPotential_',lvl(4:length(lvl)),'(lvl);']);	%calls Korek potential





R  =[4.5:.1:50];	[V,Vmorse,Vlngrng]=LiRbStates_Pot(R,lvl);						%Morse and long range potential fitted



subplot(221),plot(Rp,Vp*cnst,'.r',R,Vmorse*cnst,'g',R,Vlngrng*cnst,'b',R,V*cnst,'k')

set(gca,'xLim',[0 30],'yLim',[min(Vp) -.1*min(Vp)]*cnst,'FontSize',14,'LineWidth',1.2);

xlabel('\fontsize{14}{\itR}   ({\ita}_B)');

ylabel('\fontsize{14}{\itV}({\itR})   (cm^{-1})');

pause(.1);





[R,V,E,D]=LiRbStates_Grid(lvl);



subplot(222),plot(R,V,'b',R,3e-3*real(D).^2+ones(size(R'))*E,'g')

set(gca,'xLim',[min(R) max(R)],'yLim',[min(V) -.1*min(V)],'FontSize',14)

xlabel('\fontsize{14}{\itR}   (a_B)')

ylabel('\fontsize{14}{\itV}   (a.u.)')

% subplot(234),plot(E(2:length(E))*cnst*100*c0*1e-9/2/pi,'.r')

% set(gca,'FontSize',14)

% xlabel('\fontsize{14}{\itN}')

% ylabel('\fontsize{14}{\itE}   (GHz)')





%potential characterization



N=length(E);

for k=1:N

	Rt(k)=LiRbStates_Trn(E(k),lvl);													%outer turning point

	Bv(k)=LiRbStates_Rot(E(k),lvl);													%rotational progression simple

	Bw(k)=LiRbStates_RotW(R,D(:,k)');													%rotational progression better

	Rc(k)=sum(D(:,k).*R'.*D(:,k))/sum(D(:,k).*D(:,k));								%center-of-mass

end

subplot(234),plot(1:N,Rt,'.r',1:N,Rc,'.g')

set(gca,'FontSize',14)

xlabel('\fontsize{14}{\itN}')

ylabel('\fontsize{14}{\itR_t}   (a_B)')

subplot(235),plot(1:N,Bv*1e-6,'.r',1:N,Bw*1e-6,'.g')

set(gca,'FontSize',14)

xlabel('\fontsize{14}{\itN}')

ylabel('\fontsize{14}{\itB_v}   (MHz)')





%calculate Franck-Condon with vibrational state N=10



Psi0=D(:,10).^2;

for k=1:N

	Psi =D(:,k).^2;

	FCF(k)=sum(Psi.*Psi0)/sum(Psi.*Psi)

end

subplot(236),plot(FCF,'g')

set(gca,'FontSize',14)

xlabel('\fontsize{14}{\itN}')

ylabel('\fontsize{14}{\itFranck-Condon overlap with N=10}   (a.u.)')

end





function [Rp,Vp] = LiRbPotential_1Pi(lvl)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Energies versus R (in au) for the lowest 1Pi states of LiRb [Korek00].



M=[

%R(a0)	(1)1Pi 	(2)1Pi 	(3)1Pi 	(4)1Pi 	(5)1Pi 

4.500	-0.264855 	-0.252309 	-0.227197 	-0.215887 	-0.205977 

4.600 	-0.267549 	-0.255450 	-0.229500 	-0.219860 	-0.210060 

4.800 	-0.272860 	-0.261151 	-0.233960 	-0.227159 	-0.217609 

5.000 	-0.277899 	-0.266157 	-0.238283 	-0.233574 	-0.224277 

5.200 	-0.282479 	-0.270552 	-0.242360 	-0.239111 	-0.230051 

5.400 	-0.286497 	-0.274407 	-0.246160 	-0.243837 	-0.234981 

5.600 	-0.289916 	-0.277765 	-0.249609 	-0.247848 	-0.239157 

5.800 	-0.292745 	-0.280656 	-0.252667 	-0.251211 	-0.242621 

6.000 	-0.295027 	-0.283114 	-0.255356 	-0.253990 	-0.245467 

6.200 	-0.296808 	-0.285168 	-0.257642 	-0.256183 	-0.247776 

6.400 	-0.298153 	-0.286858 	-0.259739 	-0.257861 	-0.249636 

6.600 	-0.299119 	-0.288214 	-0.261499 	-0.259101 	-0.251057 

6.800 	-0.299770 	-0.289274 	-0.263053 	-0.259998 	-0.252160 

7.000 	-0.300155 	-0.290068 	-0.264303 	-0.260577 	-0.252936 

7.200 	-0.300333 	-0.290636 	-0.265295 	-0.260915 	-0.253448 

7.400 	-0.300346 	-0.290991 	-0.266068 	-0.261025 	-0.253712 

7.600 	-0.300241 	-0.291170 	-0.266719 	-0.260993 	-0.253822 

7.800 	-0.300051 	-0.291193 	-0.267202 	-0.260819 	-0.253777 

8.000 	-0.299807 	-0.291083 	-0.267577 	-0.260555 	-0.253642 

8.200 	-0.299531 	-0.290854 	-0.267804 	-0.260192 	-0.253364 

8.400 	-0.299240 	-0.290523 	-0.267897 	-0.259753 	-0.252945 

8.600 	-0.298952 	-0.290137 	-0.267997 	-0.259303 	-0.252611 

8.800 	-0.298668 	-0.289694 	-0.267981 	-0.258800 	-0.252147 

9.000 	-0.298393 	-0.289211 	-0.267884 	-0.258252 	-0.251560 

9.200 	-0.298134 	-0.288725 	-0.267785 	-0.257723 	-0.251098 

9.400 	-0.297885 	-0.288230 	-0.267622 	-0.257158 	-0.250539 

9.600 	-0.297651 	-0.287751 	-0.267439 	-0.256576 	-0.249939 

9.800 	-0.297428 	-0.287292 	-0.267252 	-0.255991 	-0.249377 

10.000 	-0.297215 	-0.286852 	-0.267052 	-0.255379 	-0.248766 

10.200 	-0.297014 	-0.286455 	-0.266840 	-0.254772 	-0.248195 

10.400 	-0.296823 	-0.286091 	-0.266628 	-0.254137 	-0.247595 

10.600 	-0.296642 	-0.285763 	-0.266420 	-0.253507 	-0.247028 

10.800 	-0.296470 	-0.285470 	-0.266218 	-0.252876 	-0.246477 

11.000 	-0.296305 	-0.285209 	-0.266021 	-0.252229 	-0.245904 

11.500 	-0.295928 	-0.284690 	-0.265582 	-0.250658 	-0.244530 

12.000 	-0.295597 	-0.284333 	-0.265234 	-0.249235 	-0.243342 

12.500 	-0.295307 	-0.284093 	-0.264947 	-0.247974 	-0.242119 

13.000 	-0.295055 	-0.283935 	-0.264716 	-0.246939 	-0.240888 

13.500 	-0.294838 	-0.283839 	-0.264532 	-0.246141 	-0.239804 

14.000 	-0.294654 	-0.283785 	-0.264380 	-0.245534 	-0.238784 

14.500 	-0.294494 	-0.283751 	-0.264255 	-0.245081 	-0.237935 

15.000 	-0.294362 	-0.283731 	-0.264149 	-0.244739 	-0.237246 

15.500 	-0.294249 	-0.283720 	-0.264060 	-0.244483 	-0.236726 

16.000 	-0.294151 	-0.283713 	-0.263982 	-0.244281 	-0.236313 

16.500 	-0.294075 	-0.283706 	-0.263915 	-0.244121 	-0.236005 

17.000 	-0.294008 	-0.283708 	-0.263856 	-0.243995 	-0.235778 

17.500 	-0.293952 	-0.283708 	-0.263805 	-0.243892 	-0.235605 

18.000 	-0.293905 	-0.283693 	-0.263761 	-0.243807 	-0.235471 

18.500 	-0.293864 	-0.283687 	-0.263722 	-0.243735 	-0.235360 

19.000 	-0.293831 	-0.283684 	-0.263689 	-0.243674 	-0.235274 

19.500 	-0.293802 	-0.283643 	-0.263660 	-0.243622 	-0.235200 

20.000 	-0.293777 	-0.283657 	-0.263634 	-0.243576 	-0.235136 

20.500 	-0.293757 	-0.283669 	-0.263590 	-0.243536 	-0.235020 

21.000 	-0.293739 	-0.283721 	-0.263572 	-0.243501 	-0.234976 

21.500 	-0.293724 	-0.283728 	-0.263556 	-0.243471 	-0.234934 

22.000 	-0.293711 	-0.283729 	-0.263544 	-0.243443 	-0.234893 

22.500 	-0.293699 	-0.283728 	-0.263521 	-0.243417 	-0.234827 

23.000 	-0.293690 	-0.283734 	-0.263540 	-0.243397 	-0.234866 

23.500 	-0.293681 	-0.283735 	-0.263531 	-0.243378 	-0.234834 

24.000 	-0.293674 	-0.283733 	-0.263518 	-0.243360 	-0.234797 

24.500 	-0.293668 	-0.283734 	-0.263516 	-0.243347 	-0.234773 

25.000 	-0.293662 	-0.283735 	-0.263506 	-0.243331 	-0.234738 

25.500 	-0.293657 	-0.283735 	-0.263504 	-0.243320 	-0.234724 

26.000 	-0.293653 	-0.283734 	-0.263498 	-0.243308 	-0.234700 

26.500 	-0.293650 	-0.283734 	-0.263495 	-0.243299 	-0.234679 

27.000 	-0.293646 	-0.283733 	-0.263490 	-0.243289 	-0.234660 

27.500 	-0.293643 	-0.283734 	-0.263487 	-0.243282 	-0.234643 

28.000 	-0.293641 	-0.283736 	-0.263485 	-0.243276 	-0.234631 

28.500 	-0.293638 	-0.283736 	-0.263446 	-0.243269 	-0.234549 

29.000 	-0.293636 	-0.283737 	-0.263456 	-0.243264 	-0.234551 

29.500 	-0.293634 	-0.283735 	-0.263449 	-0.243259 	-0.234530 

30.000 	-0.293632 	-0.283735 	-0.263447 	-0.243254 	-0.234519 

31.000 	-0.293629 	-0.283737 	-0.263449 	-0.243248 	-0.234512 

32.000 	-0.293627 	-0.283738 	-0.263448 	-0.243242 	-0.234503 

33.000 	-0.293625 	-0.283731 	-0.263432 	-0.243236 	-0.234449 

34.000 	-0.293624 	-0.283735 	-0.263433 	-0.243232 	-0.234454 

35.000 	-0.293622 	-0.283734 	-0.263428 	-0.243229 	-0.234438 

36.000 	-0.293622 	-0.283736 	-0.263429 	-0.243227 	-0.234442 

37.000 	-0.293621 	-0.283736 	-0.263428 	-0.243225 	-0.234452 

38.000 	-0.293619 	-0.283733 	-0.263429 	-0.243222 	-0.234445 

39.000 	-0.293619 	-0.283734 	-0.263428 	-0.243221 	-0.234438 

40.000 	-0.293619 	-0.283738 	-0.263430 	-0.243221 	-0.234454 

41.000 	-0.293618 	-0.283736 	-0.263424 	-0.243219 	-0.234437 

42.000 	-0.293618 	-0.283735 	-0.263426 	-0.243218 	-0.234437 

43.000 	-0.293617 	-0.283740 	-0.263434 	-0.243217 	-0.234450 

44.000 	-0.293615 	-0.283739 	-0.263461 	-0.243198 	-0.234512 

45.000 	-0.293616 	-0.283739 	-0.263463 	-0.243216 	-0.234513 

46.000 	-0.293616 	-0.283740 	-0.263463 	-0.243216 	-0.234513 

47.000 	-0.293616 	-0.283740 	-0.263463 	-0.243215 	-0.234513 

48.000 	-0.293616 	-0.283740 	-0.263463 	-0.243215 	-0.234512 

49.000 	-0.293616 	-0.283740 	-0.263463 	-0.243214 	-0.234511 

50.000 	-0.293616 	-0.283740 	-0.263463 	-0.243214 	-0.234511 

55.000 	-0.293616 	-0.283733 	-0.263421 	-0.243213 	-0.234426 

60.000 	-0.293615 	-0.283739 	-0.263462 	-0.243212 	-0.234509 

500.000 	-0.293614 	-0.283739 	-0.263462 	-0.243211 	-0.234508 

];



Rp  =M(:,1);

Vp  =M(:,str2num(lvl(2))+1)-M(length(M),str2num(lvl(2))+1);

end



function [V,Vmorse,Vlngrng] = LiRbStates_Pot(R,lvl)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LiRb potentials, adapts the potential in atomic units.



	Constant;

	cnst=e0^2/4/pi/E0/aB/2/pi/hP/100/c0;													%conversion a.u. to cm^-1



	[Rm,Dm,Bm,Rt,C6,C8,C10]=LiRbStates_Par(R,lvl);										%gets van der Waals parameters



	Vmorse=Dm/cnst*((1-exp(-(R-Rm)/Bm)).^2-1);											%Morse part of the potential in a.u.



	Vlngrng=-(C6./R.^6)-(C8./R.^8)-(C10./R.^10);											%long range part of the potential in a.u.



	F=exp(-(R/Rt).^10);																		%join potentials

	V=Vmorse.*F + Vlngrng.*(1-F);

end



function [Rm,Dm,Bm,Rt,C6,C8,C10] = LiRbStates_Par(R,lvl)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LiRb potentials, gets van der Waals parameters,

%  Rm : center of exchange potential			in aB

%  Dm : potential depth							in 1/cm

%  Bm : width of exchange potential			in aB

%  Rt : short-to-long range connection point	in aB

%  C6 : van der Waals coefficient				in aB^6 

%  Te : potential depth							in 1/cm

%  Re : location of potential minimum			in Angs

%  we :											in 1/cm

%  Be : rotational constant						in 1/cm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



if     sum(lvl(1:6)=='(1)1Si')==6 wch=1;

elseif sum(lvl(1:6)=='(2)1Si')==6 wch=2;

elseif sum(lvl(1:6)=='(3)1Si')==6 wch=3;

elseif sum(lvl(1:6)=='(1)3Si')==6 wch=4;

elseif sum(lvl(1:6)=='(2)3Si')==6 wch=5;

elseif sum(lvl(1:6)=='(3)3Si')==6 wch=6;

elseif sum(lvl(1:6)=='(1)1Pi')==6 wch=7;

elseif sum(lvl(1:6)=='(2)1Pi')==6 wch=8;

elseif sum(lvl(1:6)=='(1)3Pi')==6 wch=9;

elseif sum(lvl(1:6)=='(2)3Pi')==6 wch=10;

end



Pa=[																														%[Korek00 & Marinescu99]]

%Rm		Dm		Bm		Rt		C6			C8			C10			Te		Re		we		Be

6.4		4700	2.3		9.6		2545		2.34e5		2.61e7		0		3.43	194.0	.220		%(1)1Sigma+	X	%[Derevianko01 & 03]

7.84	7087	4.02	16.0	2.6744e4	2.8089e6	0			11639	4.138	119.6	.152		%(2)1Sigma+	A	

7.9    3490	3.0    15.8	-1.6506e4	1.8187e6	0			17438	4.257	113.0	.144		%(3)1Sigma+	C

9.65	279		2.5		20		3015		7.6e5		7e7			5678	5.126	41.2	.100		%(1)3Sigma+	a

7.64	3941	2.94	24.1	2.6744e4	2.8089e6	0			14719	4.064	129.9	.157		%(2)3Sigma+	c

7.35	370		1.0		16.6	-1.6506e4	1.8187e6	0			0		0		0		0			%(3)3Sigma+

7.35	1429	2.13	13.7	9.4312e3	1.7135e5	0			17205	3.813	140.5	.179		%(1)1Pi		B

7.8    1481	2.4		10.9	-2.5677e3	4.1208e5	0			19201	4.101	141.7	.155		%(2)1Pi		D

6.14	8426	2.63	11.7	9.4312e3	1.7135e5	0			10233	3.348	191.7	.232		%(1)3Pi		b

7.74	1383	2.30	13.6	-2.5677e3	4.1208e5	0			19484	4.115	39.4	.160		%(2)3Pi

];



Rm=Pa(wch,1);																							%in aB

Dm=Pa(wch,2);																							%in a.u.

Bm=Pa(wch,3);																							%in 1/aB

Rt=Pa(wch,4);																							%in 1/aB

C6 =Pa(wch,5);																						%in aB^6

C8 =Pa(wch,6);																						%in aB^8

	C10=Pa(wch,7);																						%in aB^10



	Po=[																														%[Marinescu97]

%C5			C6			C8			C10

	0			2.7415e4	8.6940e6	3.4006e9															%2S-6S 1,3Sigma+

	0			1.1736e5	8.6022e6	2.5892e9															%3S-5S 1,3Sigma+

	0			1.5217e4	5.8713e6	0 																	%2S-4D	1,3Sigma+

	0			1.3195e4	1.5770e6	0																	%2S-4D	1,3Pi	

	0			7.1308e3	-1.5547e5	0																	%2S-4D	1,3Delta

	0			2.6838e4	3.2636e7	0																	%3D-5S	1,3Sigma+

	0			2.4318e4	1.2976e7	0																	%3D-5S	1,3Pi	

	0			1.6758e4	6.8198e5	0																	%3D-5S	1,3Delta

	-3.8417e2	6.5618e3	-1.8430e6	0																	%2P-5P	1,3Delta

	1.5367e3	6.1659e3	1.7324e6	0																	%2P-5P	1,3Pi

	0			1.1244e4	-2.5382e6	0																	%2P-5P	1,3Pi

	0			4.0227e3	-2.3540e4	0																	%2P-5P	1,3Sigma-

	0			1.0107e4	-3.4418e6	0																	%2P-5P	1,3Sigma+

	-2.3050e3	-9.0214e4	3.9927e6	0																	%2P-5P	1,3Sigma+

	];

end



function Rt = LiRbStates_Trn(E,lvl)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LiRb potential, finds outer turning point.



	Rx=[.01:.01:100];

	[V,Vmorse,Vlngrng]=LiRbStates_Pot(Rx,lvl);

	[Vmin,Imin]=min(V);

	V =V (Rx>Rx(Imin));

	Rx=Rx(Rx>Rx(Imin));

	[Vmin,Imin]=min((E-V).^2);

	Rt=Rx(Imin);

end



function Bv = LiRbStates_Rot(E,lvl)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LiRb potential, finds rotational progression via course estimation.



Constant;

mred=uA/(1/6+1/87);																		%reduced mass



Rt=LiRbStates_Trn(E,lvl);																%outer turning point



Bv=hP/pi/mred./(Rt*aB).^2;																%rotational progression

end



function Bv = LiRbStates_RotW(R,Ps)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LiRb potential, finds rotational progression via course estimation.



	Constant;

	mred=uA/(1/6+1/87);															%reduced mass

	Rt=sqrt(sum(Ps.*R.^2.*Ps)/sum(Ps.*Ps));										%root mean square

	Bv=hP/pi/mred./(Rt*aB).^2;													%rotational progression

end



function [v,a] = LiRbStates_Vib(R,Psi)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LiRb potential, counts vibrational quantum number and determines scattering length.



%counts vibrational quantum number

	v=0;

	for j=2:length(Psi)

		if Psi(j)*Psi(j-1)<0

			v=v+1;

		end

	end



%determines scattering length 

	for j=length(Psi):-1:2

		if Psi(j)*Psi(j-1)<0

			a=R(j);

			break;

		end

	end

end



function [R,V,EE,DD] = LiRbStates_Grid(lvl)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fourier grid for solving a LiRb potential, potential is called from fit to [Korek00] and [Derevianko01].



Constant;

cnst=e0^2/4/pi/E0/aB/2/pi/hP/100/c0;													%conversion a.u. to cm^-1

mred=uA/(1/6+1/87);	



R  =[4:.05:30];																			%Morse and long range potential fitted

[V,Vmorse,Vlngrng]=LiRbStates_Pot(R,lvl);





N  =length(R);																			%Fourier-grid calculation

kl =[1:N];

[k,l]=meshgrid(kl,kl);

Hkl=pi^2/4/(mred/mE)/(max(R)-min(R))^2*(-1).^(k-l) .* (1 ./sin(pi*(k-l)/2/N+0e-24).^2 - 1 ./sin(pi*(k+l)/2/N+0e-24).^2);

for j=1:N

	Hkl(j,j)=0;

end

Hkk=pi^2/4/(mred/mE)/(max(R)-min(R))^2 * ((2*N^2+1)/3 - 1 ./sin(pi*kl/N+0e-24).^2) + V;

H  =Hkl+diag(Hkk);

[D,E]=eig(H);



v=0;

for j=N:-1:1

	if E(j,j)<0

		v=v+1;

		EE(v)=E(j,j);

		DD(:,v)=D(:,j);

	end

end

end